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Numerical  Simulation  of  a  Turbulent  Flow 
Through  a  Shock  Wave 


INTRODUCTION 

Renewed  interest  in  designing  a  hypersonic  aircraft  makes  essential  a 
reappraisal  of  the  various  fluid  mechanics  phenomena  typical  of  such  high  speed 
flow.  At  hypersonic  speeds  there  are  flow  phenomena  that  do  not  occur  at 
subsonic  speeds,  such  as  real  gas,  chemical  and  thermal  effects  and,  in  particular, 
the  effect  on  the  flow  of  strong  shock  waves.  Because  of  the  extremely  harsh 
environment  in  which  these  phenomena  occur,  it  is  difficult  to  investigate  them 
experimentally  and,  consequently,  a  much  greater  reliance  must  be  placed  on 
numerical  simulations  than  is  the  case  in  lower  speed  regimes. 

This  report  describes  the  results  of  research  conducted  over  a  two-year  period 
to  provide  a  foundation  for  an  extended  research  effort  into  the  interaction  between 
shock  waves  and  turbulent  flows  at  hypersonic  speeds  using  numerical  simulations. 
The  work  described  herein  was  conducted  at  transonic  speeds  in  order  to  make 
maximum  use  of  existing  knowledge  and  computational  methods  for  developing 
insight  to  the  shock/turbulence  interaction  and  to  aid  in  identifying  problems  in  the 
overall  approach. 

One  of  the  basic  problems  in  high  speed  flow  with  shock  waves  is  the 
shock/boundary  layer  interaction,  particularly  the  interaction  of  the  very  large 
pressure  gradient  with  the  turbulence.  The  lack  of  detailed  knowledge  of  this 
interaction  is  a  probable  cause  of  the  frequent  disagreement  between  numerical 
predictions  of  a  shock/boundary  layer  interaction  and  experimental  data  from  such 
interactions.  Even  if  the  shock/turbulence  interaction  itself  were  not  the  cause  of 
inaccurate  predictions,  the  uncertainty  of  the  nature  of  the  interaction  precludes  an 
accurate  identification  of  the  real  cause.  Consequently,  it  is  essential  to  understand 
the  nature  of  the  shock/turbulence  interaction. 
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It  is  expected  that  shock/turbulence  interactions  will  be  most  important  at 
hypersonic  speeds  where  the  shock  is  strongest.  However,  even  a  weak  shock  wave, 
typical  of  those  at  transonic  speeds,  should  have  a  significant  effect  on  the 
turbulence  because  of  the  very  large  pressure  gradient  across  the  shock.  The  first 
step  in  this  investigation  was,  therefore,  to  study  the  interaction  at  these  lower 
speeds,  where  the  problem  is  more  tractable,  prior  to  attempting  to  study  a 
hypersonic  shock/turbulence  interaction. 

The  primary  objectives  of  the  effort  described  herein  were  to  choose  a 
numerical  scheme  which  would  provide  calculations  of  sufficient  accuracy  to  give 
insight  into  shock/turbulence  interactions,  to  allow  extension  of  the  work  to  higher 
Mach  numbers  and  more  complex  flow  configurations,  and  to  define  a  flow 
configuration  which  could  be  used  to  study  the  interaction  between  turbulence  and 
shock  waves.  The  configuration  needed  to  be  physically  realistic,  but  also  simple 
enough  to  allow  the  identification  of  different  aspects  of  the  interaction.  The 
numerical  method  chosen  was  an  explicit,  second  order  accurate,  MacCormack 
scheme.  The  configuration  chosen  was  a  normal  shock  in  a  flow  which  was  uniform 
except  within  a  narrow  region  in  the  center  of  the  computational  domain. 
Calculations  were  performed  in  a  supersonic  flow  with  an  upstream  Mach  number  of 
1.15.  Turbulence  was  simulated  in  a  finite  portion  of  the  entrance  boundary  using 
pseudo-random  distributions  of  velocities,  density,  and  pressure  produced  by  a 
steindard  random  number  generator  function.  The  shock  wave  was  constrained  to 
be  held  at  a  fixed  position  in  the  flow  field  by  conditions  imposed  at  the  lateral 
boundaries.  Calculations  were  performed  for  the  turbulent  flow  with  and  without 
the  shock.  A  statistical  analysis  was  also  performed  by  calculating  samples  of 
various  statistical  quantities  eis  the  calculations  proceeded. 

The  remainder  of  this  report  will  describe  the  selection  of  the  numerical 
methods  used,  including  tests  of  various  computational  grid  configurations  and  the 
definition  of  the  flow  configuration.  The  particular  grid  and  computer  codes  used 
will  be  described,  along  with  the  numerical  methods  of  calculating  statistical  data 
from  calculated  turbulence.  Results  will  be  presented  for  a  sample  shock/turbulence 
interaction  and  for  a  separate  analysis  of  a  perturbed  shock  wave  based  on 
transonic  indicial  theory.  The  comparison  of  the  two  theories  will  be  shown  to 
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indicate  that  the  motion  of  the  shock  wave  due  to  turbulent  fluctuations  must  be 
accounted  for  in  order  to  accurately  model  the  passage  of  turbulence  through  a 
shock. 


NUMERICAL  METHOD 

The  choice  of  a  method  of  obtaining  numerical  solutions  to  the  governing  flow 
equations  was  made  through  considerations  of  the  nature  and  complexity  of  the  flow 
field  to  be  calculated.  Of  primary  concern  was  that  the  numerical  method  be  easily 
applied  to  the  problem  of  interest,  and  be  easily  modified  and  capable  of  extension 
to  calculate  fully  three-dimensional  flows.  Of  equal  importance  was  the  need  for 
accuracy,  both  in  time  and  space  in  order  to  be  able  to  resolve  rapid  small  scale 
fluctuations. 

Since  the  calculations  to  be  made  were  intended  to  provide  information  about 
the  fundamental  nature  of  the  flow,  the  need  for  accuracy  and  ease  of  use  of  the 
method  were  also  accompanied  by  the  need  for  minimal  artificial  factors  attributable 
to  the  numerical  method  itself.  Examples  of  such  artificial  factors  are  artificial 
dissipation,  shock  smearing,  that  is,  spreading  of  the  shock  over  several 
computational  points,  and  numerical  oscillations.  Also,  the  method  must  be  able  to 
calculate  the  correct  shock  speed  and  be  well  suited  to  a  time  accurate  calculation. 
Finally,  the  computational  method  should  be  efficient,  allowing  many  calculations  to 
be  made  at  minimum  cost. 

Three  types  of  numerical  methods  were  examined:  artificial  dissipation  methods, 
linear  hybridization  methods,  and  Godunov  methods.  Artificial  dissipation  methods, 
such  as  generalized  Lax-Wendrof  methods,  including  MacCormack  schemes,  were 
found  to  be  desirable  because  of  simplicity  of  application.  Such  methods  are 
especially  successful  when  the  artificial  viscous  effects  can  be  mostly  confined  to  a 
well-defined  finite  region,  such  as  the  high  gradient  region  associated  with  a  shock 
wave. 


Linear  hybridization  methods,  such  as  flux  corrected  transport  or  total 
variational  diminishing  (TVD)  methods  have  the  advantage  of  improved  resolution  of 
flow  discontinuities  when  a  uniform  grid  is  used.  However,  the  methods  require 
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double  the  work  of  artificial  dissipation  methods  and  are  clearly  limited  by  the 
resolving  power  of  the  low  order  fluxes.  Also,  there  is  difficulty  in  devising 
appropriate  weight  factors  for  the  low  order  and  high  order  fluxes. 

Godunov's  approach  avoids  artificial  dissipation  and  is  highly  accurate  when  a 
higher  order  difference  scheme  involving  explicit  nonlinearity  is  used.  However,  the 
use  of  a  Riemann  solver  introduces  a  great  deal  of  complexity  into  a  difference 
scheme.  Use  of  such  a  method  for  the  problem  of  interest  herein  would  involve  a 
major  code  development,  a  formidable  research  program  in  itself.  Such  an  approach 
was  clearly  unsuitable  for  the  present  study. 

The  computational  method  chosen  for  the  present  study  was  the  explicit 
MacCormack  scheme.  The  choice  was  based  on  three  favorable  characteristics  of 
the  method:  (l)  Simplicity;  a  two-dimensional  scheme  could  be  developed  for 
preliminary  work  and  could  then  be  easily  extended  to  three  dimensions;  (2) 
Relatively  good  resolution  of  unsteady  shocks;  (3)  An  estimate  of  implicit  artificial 
dissipation  is  available,  allowing  an  evaluation  of  the  quality  of  the  numerical 
results. 

A  two-dimensional  Navier-Stokes  solver  was  developed  using  the  standard 
MacCormack  scheme  algorithm  in  a  clustered  cartesian  coordinate  system.  An  H- 
grid  was  chosen  to  avoid  coordinate  transformations  for  comparison  of  Reynolds 
stresses  before  and  after  shocks.  The  shock  capturing  performance  of  the  code  was 
tested  by  application  to  a  flow  field  with  an  oblique  shock.  The  results  were 
compared  to  results  obtained  from  a  similar  code  developed  by  J.  Shang  (Ref.  1). 
Examples  of  the  comparison  of  the  two  codes  axe  shown  in  Figure  1.  In  Figure  la 
and  b  plots  of  constant  density  contours  for  an  oblique  shock  at  freestream  Mach 
number  of  2.0  are  shown.  The  code  of  Reference  1  is  designated  Code  1  in  the 
figure,  and  the  code  developed  in  the  present  work  is  Code  2.  Both  codes  produce 
sbnilax  solutions  from  the  perspective  given  by  the  density  contours.  In  Figure  Ic 
and  d,  the  density  along  a  line  through  the  shock  wave  is  shown  for  the  two  codes. 
Clearly,  both  codes  suffer  somewhat  from  the  coarseness  of  the  grid  and  numerical 
dispersion  effects.  However,  the  new  code  produces  a  solution  with  reduced 
oscillations.  The  effects  on  the  code's  shock  capturing  performance  of  grid 
refinement,  kinematic  viscosity  and  second  and  fourth  order  artificial  dissipation  were 
tested,  and  results  will  be  presented  subsequently. 
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FLOW  CONFIGURATION 


In  choosing  a  configuration  for  calculations  of  shock  wave  turbulence 
interactions,  it  was  decided  initially  that  a  flow  field  was  required  that  would  allow 
observation  of  the  direct  effect  of  the  shock  wave  on  the  turbulence  and,  conversely, 
of  the  turbulence  on  the  shock.  The  flow  field  should  be  clearly  defined  and  allow 
the  researcher  as  much  control  as  possible  over  the  essential  elements  of  the  flow 
through  manipulation  of  the  inflow,  outflow,  and  boundary  conditions.  Thus,  it  was 
considered  important  to  avoid  such  things  as:  curved  geometry,  which  would 
require  a  curvilinear  coordinate  system;  oblique  shocks,  which  would  possibly 
undergo  bifurcation,  thus  obscuring  the  effects  of  interest;  shock  reflections  from 
boundaries,  which  would  create  a  complex  flow  field,  thus  making  more  difficult  the 
interpretation  of  shock-turbulence  interaction  effects;  and  no-slip  boundary  conditions 
which  would  create  a  shear  flow  and  lead  to  development  of  a  more  complex  field 
of  turbulence. 

The  first  configuration  chosen  for  examination  was  a  flow  in  a  duct  with  a 
bump  on  one  wall.  The  incoming  flow  was  subsonic,  and  a  transpiration  condition 
was  to  be  applied  to  the  lower  wall  to  generate  a  normal  shock,  with  a  symmetry 
condition  at  the  upper  boundary.  The  exit  flow  would  also  be  subsonic. 
Unfortunately,  such  a  configuration  violates  several  of  the  criteria  mentioned 
previously  as  to  be  avoided.  For  example,  the  boundary  condition  on  the 
streamwise  velocity  component  on  the  lower  boundary  is  a  no-slip  condition.  In 
addition  to  the  creation  of  shear,  such  a  condition  leads  to  slow  numerical 
convergence  because  of  the  grid  refinement  required  for  proper  resolution  of  the 
near-wall  flow.  The  no-slip  condition  was  necessary  because  the  equations  being 
solved  were  the  full  Navier-Stokes  equations.  If  the  equations  could  have  been 
reduced  to  the  Euler  equations  by  neglecting  the  viscous  terms,  the  no-slip  condition 
would  not  have  been  necessary.  However,  it  was  considered  important  to  include 
the  viscous  terms  since  the  flow  was  expected  to  be  subject  to  viscous  dissipation. 
Another  fault  of  the  duct  flow  configuration  was  that  the  shock  wave  location  was 
not  precisely  controllable,  nor  was  the  shock  perfectly  normal,  even  without 
turbulence.  Finally,  the  shock /turbulence  interaction  region  was  not  clearly  defined. 
All  these  considerations  led  to  rejecting  the  duct-with-a-bump  as  the  configuration  of 
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choice,  at  least  for  the  preliminary  study.  In  retrospect,  some  of  the  objectionable 
conditions  were  found  to  exist  to  some  degree  in  the  configuration  that  was  actually 
used,  so  that  the  duct-with-a-bump,  or  flow  over  an  airfoil,  remains  a  possibility  for 
further  study. 

The  configuration  which  was  used  was  a  simpler  one  in  which  the  walls  of  the 
duct  were  removed,  and  a  normal  shock  wave  was  created  in  a  uniform  flow  field 
which  contained  turbulence  in  a  central  region.  The  incoming  flow  was  slightly 
supersonic  (M  ==  1.15).  Inviscid  boundary  conditions  were  used  at  the  top  and 
bottom  boundaries.  The  initial  flow  field  was  specified  as  a  uniform  flow  with  a 
normal  shock,  thus  modeling  a  flow  in  a  shock  tube  in  which  the  observer  is 
moving  with  the  shock.  The  location  of  the  shock  was  fixed  in  space  by  keeping 
the  solution  at  the  top  and  bottom  boundaries  constant  and  specifying  the  pressure 
at  the  downstream  boundary  to  be  the  theoretical  pressure  for  a  normal  shock  from 
one-dimensional  Rankine-Hugoniot  theory.  Some  refinement  of  the  solution  was 
achieved  by  running  a  steady  flow  solution  in  which  the  incoming  flow  was 
constant,  and  the  top  and  bottom  boundary  conditions  were  periodically  revised  by 
setting  the  flow  variables  equal  to  their  values  at  the  centerline.  In  this  way,  after 
several  adjustments,  a  steady  flow  field  was  achieved  which  was  the  same 
everywhere,  including  the  boundaries,  with  a  shock  jump  which  spread  over  several 
grid  points  in  the  streamwise  direction,  the  number  of  grid  points  depending  upon 
the  magnitude  of  the  real  and  artificial  viscosities. 

The  generation  of  turbulence  in  the  flow  field  will  be  discussed  in  a 
subsequent  section.  It  is  noted  that  the  turbulence  was  specified  in  a  central 
section  of  the  inflow  boundary.  Disturbances  from  the  turbulent  flow  field 
propagated  along  Mach  lines  toward  the  top  and  bottom  boundaries  of  the 
configuration.  The  width  of  the  computational  space  was  such  that  these 
disturbances  intersected  the  shock  wave  before  they  intersected  the  boundaries. 
Thus,  the  boundary  conditions  should  represent  realistically  the  steady  propagation 
of  a  shock  wave  in  a  uniform  flow  field  encountering  a  finite  region  of  turbulence 
with  the  observer  moving  along  with  the  undisturbed  shock  wave. 

The  development  of  the  flow  field  in  the  steady  case  can  be  seen  in  Figures  2 
and  3.  In  Figure  2,  the  convergence  history  of  the  maximum  change  in  the 
numerical  solution  is  shown.  When  the  boundary  conditions  change,  every  200  time 
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steps,  there  is  a  jump  in  the  residual,  which  then  proceeds  to  approach  an 
asymptote.  The  asymptotic  level  is  reduced  each  time  the  boundary  conditions  are 
refined.  After  four  such  refinements,  as  shown  in  the  figure,  the  residual  remained 
constant.  It  is  noted  that  the  particular  numerical  scheme  being  used  in  these 
calculations  was  an  explicit  MacCormack  scheme,  for  which  the  asymptotic  value  of 
the  residual  is  a  function  of  the  time  step.  Use  of  an  implicit  formulation  would 
eliminate  this  effect.  In  Figure  3,  the  profile  of  the  Mach  number  variation 
through  the  shock  wave  is  shown.  The  effects  of  second  and  fourth  order  artificial 
dissipation  and  of  Reynolds  number  are  indicated.  The  first  case  is  a  Reynolds 
number  of  10,000  and  no  artificial  dissipation.  The  results  shown  in  Figure  3a 
indicate  that  the  shock  is  spread  over  about  6  grid  points,  and  the  numerical 
scheme  produces  some  oscillations  both  upstream  and  downstream  of  the  shock 
jump.  In  Figure  3b  some  second  order  dissipation  has  been  added,  producing  no 
discernible  change  in  the  solution.  Adding  some  fourth  order  dissipation  alone,  on 
the  other  hand,  reduces  substantially  the  oscillations  downstream  of  the  shock  while 
leaving  unaffected,  or  possibly  slightly  increasing,  the  upstreaun  oscillations  (Fig.  3c). 
Including  both  second  and  fourth  order  dissipation  apparently  is  no  better  than  the 
second  order  alone,  as  indicated  in  Figure  3d.  Reynolds  number  effects  are  quite 
large  as  shown  in  Figure  3e  and  f,  where  the  Reynolds  number  varies  by  an  order 
of  magnitude  above  and  below  that  of  Figure  3d.  The  higher  Reynolds  number 
(lower  viscosity)  clearly  reduces  the  spreading  of  the  shock  over  the  grid,  but 
increases  the  amplitude  of  the  oscillations.  The  lower  Reynolds  number  increases 
the  spreading  and  reduces  the  oscillation  amplitude. 

Clearly,  the  choice  of  the  parameter  values  which  govern  the  quality  of  the 
numerical  solution  is  the  set  of  values  which  give  the  best  resolution  of  the  shock 
jump  with  the  least  oscillations  for  a  given  Reynolds  number.  The  choice  of 
Reynolds  number  is  somewhat  arbitrary.  However,  considerations  of  the  desire  to 
resolve  a  wide  range  of  scales  of  turbulence  lead  to  a  value  which  allows  the  shock 
to  be  captured  in  as  few  points  as  possible  commensurate  with  the  minimization  of 
the  oscillations.  With  these  considerations,  both  second  and  fourth  order 
dissipation  parameters  were  used  for  subsequent  calculations,  with  Re  =  10,000. 
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GRID  DEFINITION 


The  computational  grid  used  in  preliminary  calculations  was  uniform  in  the 
lateral  directions  (parallel  to  the  shock  wave)  and  exponentially  stretched  in  the 
streamwise  direction  downstream  of  the  shock.  While  this  grid  configuration  allowed 
the  efficient  development  of  the  computational  techniques  for  producing  a  turbulent 
flow  and  computing  statistics,  accuracy  in  calculating  the  turbulence  required  a 
uniform  grid  since  that  would  allow  the  minimum  resolvable  scales  to  be  the  same 
over  the  entire  computational  space.  A  uniform  grid  was  therefore  used  in  all 
subsequent  calculations. 

The  computational  space  is  shown  in  Figure  4.  For  most  of  the  calculations, 
the  grid  dimensions  were  0.1  x  0.1  with  each  axis  divided  into  80  segments.  Some 
calculations  were  made  with  a  160x160  point  grid  to  test  the  effect  of  the  grid. 
The  x  axis  defines  the  flow  direction,  the  y  axis  defines  the  lateral  direction.  The 
flow  enters  the  computational  space  from  the  left,  at  x  =  0.  The  region  where 
turbulence  is  introduced  is  along  the  y  axis  in  the  center  of  the  inflow  boundary. 

STATISTICAL  ANALYSIS 

Understanding  a  turbulent  flow  requires  examination  of  several  kinds  of  data. 
Instantaneous  as  well  as  statistical  information  can  provide  useful  insights. 
Statistical  quantities  can  include  a  wide  variety  of  data,  from  simple  averages  of 
quantities  measured  (or  calculated)  at  a  single  point  in  space  or  time  to  correlations 
of  quantities  measured  at  several  points  and  combinations  of  space  and  time.  In 
this  investigation,  the  statistical  analysis  of  the  flow  dealt  primarily  with 
instanteineous  velocity,  pressure  and  density  distributions,  and  statistical  averages  of 
velocity,  pressure  and  density  fluctuations,  and  products  of  fluctuations  defined  at  a 
point.  Some  two-point  correlations  in  both  time  and  space  were  also  calculated. 

Basic  statistical  quantities  were  calculated  by  averaging  over  a  specified 
ensemble  of  grid  points  and  time  steps.  The  main  requirement  of  the  statistical 
analysis  is  that  the  quantities  included  in  the  ensemble  be  independent.  The 
turbulence  under  study  was  produced  by  an  input  distribution  of  quantities 
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randomly  distributed  in  time  and  along  the  y  axis  at  the  inflow  boundary.  The 
solution  generated  by  this  inflow  was  assumed  to  be  homogeneous  with  respect  to 
variations  in  the  y  direction  within  a  range  defined  by  the  specified  inflow.  The 
basic  statistical  ensemble  at  each  x  station  was  defined  to  consist  of  data  at  grid 
points  along  lines  parallel  to  the  y  axis  for  all  time  steps.  In  order  to  avoid  edge 
effects,  of  N  grid  points  included  in  the  inflow  turbulence  specification,  N-4  points 
were  included  in  the  statistical  ensemble.  Statistical  quantities  were  defined  by 


1  L  T 

<f>  =  /  J  f(x,t)dxdt 

LT  0  0 


or,  for  a  finite  numerical  ensemble, 


<f> 


1  J-1  N 

-  E  E  f. 

N(J-2)  j=2  n=l 


(1) 


(2) 


where  f  is  a  quantity  of  interest,  i.e.,  u,  v,  p,  p,  uv,  pu,  ^u,  etc.,  J  is  the  number 
of  y-axis  points  at  which  the  turbulence  is  introduced  at  the  inflow  boundary,  and 
N  is  the  number  of  time  steps  included  in  the  ensemble.  Values  of  all  four  basic 
flow  quantities,  u,  v,  p,  p  were  accumulated  at  36  y  values  between  y  =  0.025  and 
0.075  in  order  to  maximize  the  amount  of  data. 


REPRESENTATION  OF  TURBULENCE 

Turbulence  is  defined  as  an  irregular  condition  of  flow  in  which  the  various 
quantities  show  a  random  variation  with  time  and  space  coordinates,  so  that 
statistically  distinct  average  values  can  be  discerned  (Ref.  2).  Frequently  "pseudo¬ 
turbulence"  is  used,  referring  to  a  flow  field  with  a  regular  pattern  that  shows  a 
distinct  constant  periodicity  in  time  and  space,  but  produces  certain  statistical 
properties  of  a  turbulent  flow.  In  recent  numerical  simulation  studies  of  turbulence, 
a  typical  technique  is  to  begin  a  time-dependent  solution  of  the  Navier-Stokes 
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equations  with  a  flow  field  with  fluctuations  of  the  velocities  produced  by  various 
kinds  of  "pseudo-random"  number  generators  (Ref.  3).  After  some  integration  time, 
it  is  found  that  the  solution  has  evolved  into  a  state  that  resembles  an  actual 
turbulent  flow. 

In  the  present  study,  the  realistic  simulation  of  turbulence  was  of  less 
importance  than  the  accuracy  of  the  computational  method  and  the  definition  of  a 
flow  configuration  that  would  allow  shock/turbulence  interaction  effects  to  be 
discerned.  In  the  first  place,  the  preliminary  flow  configuration  was  two- 
dimensional,  thus  precluding  realistic  simulation  of  a  three-dimensional  field  of 
turbulent  fluctuations.  Secondly,  it  was  not  desirable  to  devote  effort  to  perfecting 
the  turbulence  simulation  until  the  fundamental  value  of  the  overall  approach  could 
be  demonstrated. 

As  a  preliminary  method  of  simulating  a  turbulent  flow  field,  the  two  velocity 
components,  and  the  pressure  and  density  were  prescribed  at  each  point  of  a  range 
of  y  values  at  the  inflow  boundary  using  a  series  of  numbers  produced  by  a 
standard  random  number  generator  which  produced  pseudo-random  values  with  a 
Gaussian  distribution  between  0  and  1.  These  numbers  were  then  modified  to 
produce  some  correlation  between  the  values  at  successive  time  steps  by  the 
expression 

i 

S  =  CS  1  +  (1  -  C^)^  R  (3) 

n  n-1  ^  n  ' 

where  Rn  is  the  sequence  of  independent  random  numbers,  Sj,  is  the  next  value  of 
the  correlated  sequence,  Sn.x  is  the  previous  value,  and  C  is  the  correlation 
coefficient.  In  this  work,  the  correlation  coefficient  wais  taken  to  be 


C  =  0.98 


(4) 


This  value  was  somewhat  arbitrarily  chosen  so  that  the  autocorrelation  of  the  Sn 
sequence  decre2ises  to  0.5  in  about  30  time  steps.  These  quantities  were  then  used 
to  define  fluctuations  of  the  flow  field  by  the  relations 
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(5) 


P  =  +  aS^) 


u  =  u  (1  +  bS^) 

(6) 

o3 

(7) 

P  =  Poo(l  ^ 

(8) 

where  the  superscript  on  the  random  sequence  identifies  a  different  sequence  for 
each  variable,  each  sequence  independent  from  the  others.  It  is  recalled  that  these 
expressions  are  applied  as  a  function  of  time  at  each  point  in  a  range  of  y  values 
at  the  inflow  boundary.  Thus,  while  the  time  variation  of  the  fluctuations  at  each 
y  location  is  not  completely  random  because  of  the  correlation  coefficient,  C,  the  y 
variation  is  more  random  due  to  the  relative  independence  of  the  terms  of  the 
sequence.  A  simple  three  point  smoothing  was  applied  to  the  inflow  conditions  in 
the  y-direction  to  help  eliminate  numerical  problems  created  by  large  changes 
between  the  specified  values  of  the  flow  variables  at  adjacent  y  locations  in  the 
incoming  flow. 

RESULTS  FOR  NORMAL  SHOCK/TURBULENCE  INTERACTIONS 

Calculations  were  performed  for  three  cases.  First,  the  solution  for  a  uniform 
steady  flow  with  a  shock  wave  was  calculated.  This  provided  the  initial  conditions 
for  the  next  case,  a  flow  with  turbulence  specified  as  described  previously  using 
random  numbers  in  a  finite  region,  forming  a  "jet"  of  turbulent  flow  in  the  center 
of  the  uniform  flow.  The  side  boundary  conditions  for  that  flow  were  the 
conditions  at  the  final  time  step  of  the  non-turbulent  flow.  Thus,  the  converged 
steady  flow  solution  provided  both  the  initial  conditions  and  the  boundary 
conditions  for  the  subsequent  turbulent  solution.  It  also  provides  an  estimate  of  the 
error  associated  with  the  numerical  procedure.  The  third  case  calculated  was  a  flow 
like  the  second  case  except  that  the  flow  was  uniform  everywhere  initially,  with 
supersonic  inflow  as  in  both  previous  cases,  but  with  no  shock  jump  on  the  side 
boundaries  and  with  all  independent  variables,  including  the  downstream  pressure. 
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calculated  by  extrapolation  at  the  downstream  boundary.  The  second  and  third 
cases  thus  provided  a  comparison  betw'een  a  turbulent  flov/  in  a  supersonic  slreniu 
and  a  turbulent  flow  passing  through  a  shock  wave. 


The  calculated  variables  were  normalized  with  respect  to  free  stream  quantities, 
with  the  resulting  steady  state  conditions  at  the  inflow  boundary  as  follows; 


Uoo 

=  1.15 

Voo 

=  0.0 

P<n 

=  1.0 

Poo 

=  1/7 

7 

=  1.4 

The  Mach  number  distribution  for  the  steady  solution  has  already  been  shown 
in  Figure  3d.  The  corresponding  pressure  distribution  is  shown  in  Figure  5. 
Typical  inflow  perturbations  are  shown  in  Figure  6.  In  Figure  6a,  the  p  and  pn 
fluctuations  as  a  function  of  time  at  the  centerpoint  of  the  inflow  boundary  are 
presented,  while  in  Figure  6b  the  variation  with  y  at  a  typical  time  step  is  shown. 
An  indication  of  the  effect  of  viscosity  on  the  turbulence  is  shown  in  Figure  7 
where  the  y  variation  of  p  and  pxi  at  an  x-station  near  the  shock  wave  is  seen  to 
be  considerably  smoother  thrm  at  the  inflow  boundary.  The  primary  features  of  the 
previous  plot  are  still  discernible  at  this  station,  but  the  smaller  features  are 
smoothed  a  great  deal. 


For  the  analysis  of  turbulence,  the  statistical  averages  of  fluctuating  quantities 
are  obtained  by  subtracting  the  mean  values.  For  any  collection  of  values  of  a 
function,  the  mean  value  is  defined  by  Equation  (l)  or  Equation  (2).  By  definition, 
then,  the  average  of  the  fluctuations  about  the  mean  value  is  zero.  That  is. 


<£'>  =  <£-  <£»  =  0 


(9) 
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Thus,  products  of  fluctuating  quantities  can  be  calculated  eis  follows 


<f'g'>  -  <(f  -  <f»  (g  -  <g>)>  =  <fg>  -  <f><g>  (10) 

In  these  definitions,  it  is  important  to  note  that  the  mean  value  defined  by 
Equation  (2)  depends  upon  the  number  of  values  used  in  the  calculation.  Also,  the 
presence  of  serial  correlation,  or  dependence  between  successive  terms  of  a  random 
sequence,  can  produce  a  bias  in  a  finite  sequence.  Given  a  mean  flow  and  random 
fluctuations  about  the  mean,  it  is  possible  to  achieve  an  average  value  for  the  mean 
flow-plus-fluctuations  that  equals  the  given  mean  flow  within  a  finite  uncertainty. 
The  uncertainty  can  be  reduced  to  a  small  value,  but  can  never  be  zero.  In  test 
calculations  with  random  numbers,  it  was  concluded  that  at  leaist  10^  samples  were 
needed  to  achieve  an  acceptable  level  of  uncertainty.  For  example.  Figure  8  shows 
the  mean  velocity  fluctuation  as  given  by  Equation  (2)  with  the  individual  values 
given  by  the  random  function  part  of  Equation  (6).  The  mean  is  presented  as  a 
function  of  the  nximber  of  values  in  the  ensemble  for  up  to  10®  points.  The  figure 
shows  that  for  ensembles  containing  less  than  50,000  points,  the  mean  value  is 
strongly  dependent  upon  the  size  of  the  ensemble.  For  larger  ensembles,  the  mean 
value  approaches  zero,  but  still  can  fluctuate  considerably.  For  the  calculations 
performed  on  the  80x80  grid  in  this  work,  values  of  the  flow  quantities  were  taken 
from  36  y  locations  at  each  time  step.  To  assure  a  reasonable  accuracy  in  the 
statistics,  a  total  of  8000  time  steps  were  calculated,  making  288,000  values  available 
for  analysis.  Test  calculations  for  2000  and  4000  time  steps  produced  results  in 
close  agreement  with  the  longer  runs. 

Some  of  the  statistics  of  the  calculated  turbulence  are  shown  in  Figure  9-12. 
In  Figure  9,  the  turbulence  intensities,  i.e.,  the  rms  values  of  the  fluctuations  in  the 
u  and  V  velocity  components,  are  shown.  The  calculated  distributions  of  the 
quantities  for  the  cases  with  and  without  a  shock  wave  are  superimposed. 
Upstream  of  the  shock,  both  solutions  are  identical,  whereas  downstream  of  the 
shock,  there  are  measurable  differences.  Both  the  streamwise  (normal  to  the  shock) 
and  the  spanwise  (parallel  to  the  shock)  turbulence  intensity  components  appear  to 
be  slightly  increased  upon  passage  through  the  shock.  The  large  spike  in  the 
intensities  at  the  location  of  the  shock  is  related  to  the  shock  capture  region.  In 
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general,  the  shock  region  is  a  nonconservative  region  in  which  the  numerical  solution 
is  meaningless,  and,  accordingly,  no  special  significance  should  be  attached  to  the 
shock  capture  region,  but  only  to  the  beginning  and  end  points  of  the  region  (Ref. 
4).  This  suggests  that  the  shock  capture  region  in  Figures  9-12  can  be  ignored  for 
the  present  analysis.  The  extent  of  the  shock  capture  region  is  best  determined  by 
examining  the  variation  of  the  mean  flow  through  the  shock.  The  distribution  of 
the  mean  pressure,  shown  in  Figure  13,  indicates  that  the  shock  region  consists  of 
seven  grid  points,  ranging  from  x  =  0.0275  to  x  =  0.0350.  The  turbulence 
intensities  shown  in  Figure  9  undergo  large  fluctuations  in  the  shock  region  and  do 
not  return  to  the  shock-free  variation  after  the  shock.  The  increase  in  the 
turbulence  intensity  downstream  of  the  shock  is  considered  to  be  due  to  the 
interaction  between  the  shock  wave  and  the  turbulence.  Figure  10  shows  that  there 
is  also  a  density  effect  and  a  small  effect  on  the  pressure  fluctuations,  while  Figure 
11  indicates  that  the  turbulence  kinetic  energy  is  increeised  by  the  shock  wave  and 
Figure  12  shows  an  increase  in  the  product  <^'u'>  from  slightly  negative  ahead  of 
the  shock  wave  to  0.0003  downstream  of  the  shock,  with  a  slightly  smaller 
magnitude  effect  on  the  product  <p'u'>. 

In  view  of  the  smallness  of  the  numerical  values  of  the  changes  in  the  density- 
velocity  products  and  the  turbulence  kinetic  energy,  a  note  about  the  accuracy  of 
the  numerical  calculations  is  in  order.  The  algorithm  used  for  the  calculations  is 
formally  of  second  order  accuracy.  This  means  that  the  solution  error  is  of  the 
order  of  (Ax)  ,  or,  for  the  grid  used  in  the  calculations  shown  in  Figures  9-13,  the 
order  of  10  .  Another  way  to  assess  the  accuracy  of  the  calculations  is  to  examine 
the  solution  for  the  case  where  no  turbulence  was  present.  In  this  example,  the 
exact  solution  is  known  and  the  effect  can  be  easily  estimated  by  comparison  of  the 
computed  and  theoretical  results.  The  variation  of  the  numerical  solution  for  the 
pressure  variation  along  the  centerline  of  the  computational  region  was  shown  in 
Figure  5.  The  oscillations  shown  in  that  figure  are  due  to  the  numerical  scheme, 
and  are  the  response  to  the  shock  discontinuity.  Another  source  of  error  of  interest 
in  this  analysis  is  the  random  error  associated  with  the  numerical  calculations. 
This  error  can  be  quantified  by  applying  the  statistical  calculations  to  +he  steady 
flow  case.  The  error  is  indicated  by  the  size  of  the  "turbulent  intensity"  or  rms 
fluctuation  of  the  solution  about  the  mean.  The  results  for  the  rms  velocity 
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components  and  the  rms  pressure  and  density  fluctuations  are  presented  in  Table  I. 

The  results  indicate  that  the  maximum  error  outside  the  shock  capture  region  is 
-5 

approximately  5x10’  .  Recall  that  the  rms  fluctuations  as  indicated  in  Figures  9 

“2 

and  10  are  of  the  order  of  5x10’  ,  three  orders  of  magnitude  larger.  The  variation 
of  the  mean  velocities,  density,  pressure,  and  mass  flux  for  the  steady  flow  case  is 
listed  in  Table  II,  along  with  the  theoretical  values  for  a  normal  shock  wave.  The 
maximum  errors  in  the  mean  quantities  occur  ^ear  the  shock  due  to  the  oscillations 
of  the  numerical  scheme. 


Consider  now  the  velocity  and  density  product.  If  each  quantity  is  taken  to 
be  composed  of  an  exact  mean  value,  an  error,  and  a  turbulent  fluctuation,  that  is 

u  =  Uq  +  Ug  +  u'  +  Ug'  (11) 

P  =  Po  ^  ^  P'  ^  Pe‘  (^2) 


where  (  )o  =  the  exact  solution 


(  )e  =  the  error 


and 


(  ) '  =  the  turbulent  fluctuations,  then  the  mean  values  over  the  entire 

ensemble  of  values  calculated  are 


<u>  =  +  <u^> 

o  e 


(13) 


and 


<P>  =  Po  *  <Pe> 


(14) 


and 


<^u>  -  </j><u>  =  </j'u'>  +  Error 


(15) 


where 


Error  = 


(16) 
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The  largest  error  in  Ug'  is  due  to  the  oscillations  from  the  shock  capture. 
These  oscillations  are  spatial  and  are  relative  to  the  shock.  If  the  shock  does  not 
move,  then  there  is  no  contribution  to  Ug‘.  If  the  shock  does  move,  then  the 
maximum  error  can  be  estimated  by  assuming  the  error,  E,  is  proportional  to  the 
magnitude  of  the  oscillations,  Eq,  and  is  given  by 

E  =  2E  dX  /AX 
o  s' 

where  fiXg  is  the  shock  movement  and  fiX  is  a  grid  cell;  the  wavelength  of  the 
oscillation  is  2AX.  Taking  the  rms.  value  of  dXg  as  a  guide,  and  noting  that 

E  =:  2(10-^),(dX^)^_^^  =  8(10-2) 

it  is  found  that 

E  =  10"^ 
or 


Hence  the  overall  error  in  <^'u'>  is  of  the  order  of  1%. 

Another  perspective  on  the  errors  in  the  calculations  is  afforded  by  the 
quantity  </Ju>.  For  an  exact  one-dimensional  steady  shock  solution,  this  quantity 
should  be  constant.  The  results  in  Table  II  for  the  case  with  no  turbulence 
indicate  that  the  oscillations  of  the  numerical  solution  can  correspond  to  a  jump  in 
<fiu>  of  the  order  of  *3.5  x  10  depending  on  where  the  data  are  taken. 
Comparing  Figures  14  and  15  showing  the  correlation  <pu>  for  the  cases  without 
and  with  turbulence  respectively,  it  appears  that  the  error  in  that  quantity  is  nearly 
the  same,  independent  of  the  presence  of  turbulence.  The  situation  is  clarified  by 
comparing  Figure  14  with  Figure  16  which  shows  the  difference  between  the  <pu> 
distribution  with  turbulence  and  a  shock  (Fig.  15)  and  without  a  shock  as  in 
Figure  17.  The  fact  that  the  oscillations  are  of  the  same  order  as  the  observed 
jump  in  the  turbulent  <p'u'>  raises  the  question  of  whether  the  turbulence  is 
affected  by  the  oscillations.  Examining  the  solutions  discussed  previously  with 
reference  to  Figures  5  and  9-13,  the  oscillations  present  in  the  instantaneous 
pressure  distribution  (Fig.  5)  are  much  less  obvious  in  the  averaged  turbulent  flow 
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(Fig.  13).  Also,  while  still  present,  the  numerical  oscillations  in  <p'u‘> 
approaching  the  shock  (Fig.  12)  are  much  smaller  in  amplitude  than  the  change  in 
</3'u'>  across  the  shock,  and  the  smoothness  of  indicates  that  it  is  not 

affected  significantly  by  the  oscillatory  numerical  error  seen  in  the  </5u>  and 
<p><\i>  distributions  (Figs.  15  and  18).  Thus,  it  is  concluded  that  the  small 
changes  found  in  the  turbulence/shock-wave  interaction  are  due  to  flow  interactions 
and  not  to  numerical  error. 

Further  confirmation  of  this  conclusion  was  obtained  by  calculating  the  flow  on 
a  finer  grid  with  160  points  on  each  axis.  A  calculation  of  3300  time  steps  with 
the  turbulence  defined  at  81  points  in  the  y  direction  provided  the  results  shown  in 
Figures  19-22.  The  shock  region  as  defined  by  the  transition  from  upstream 
conditions  to  downstream  conditions  and  exemplified  by  the  mean  pressure 
distribution  in  Figure  19  is  slightly  narrower  than  for  the  coarser  grid.  The  shock 
extends  from  x  =  0.0275  to  0.03375,  over  11  grid  points.  Quantitatively,  the 
solutions  on  the  two  grids  cannot  be  directly  compared,  since  the  turbulence  inflow 
conditions  were  defined  at  each  point  of  the  inflow  boundary  in  both  cases  and  are 
therefore  different  in  the  two  cases.  Also,  the  time  step  used  for  the  fine  grid 
calculations  was  half  that  used  for  the  coarser  grid.  Thus,  the  turbulence  generated 
on  the  fine  grid  was  of  a  finer  scale,  both  spatially  and  temporally  than  that  on 
the  coarser  grid.  However,  the  solutions  for  both  grids  lead  to  the  same  conclusion, 
that  the  turbulence  is  affected  by  the  shock  wave.  The  turbulence  intensities  shown 
in  Figure  20  are  of  the  same  order  as  those  in  Figure  9.  The  change  in  the  v' 
component  downstream  of  the  shock  is  nearly  the  same  in  both  cases.  The  u' 
component,  on  the  other  hand,  appears  to  undergo  slightly  less  change  for  the  fine 
grid.  The  density  and  pressure  effects  also  are  similar  to  those  on  the  previous 
grid.  Comparing  Figures  21  and  10,  the  pressure  effect  is  nearly  the  same  on  both 
grids,  while  the  density  effect  is  somewhat  less  on  the  fine  grid.  Turbulence  kinetic 
energy.  Figure  22,  is  slightly  decreased  on  the  fine  grid  from  that  shown  in  Figure 
11.  However,  the  basic  effect  of  an  increase  over  the  no-shock  solution  is  present 
in  both  solutions.  Finally,  the  fine  grid  produces  a  jump  in  both  the  mean  p'u' 
correlation  and  the  mean  ^'u'  correlation  that  is  smaller  on  the  fine  grid, 
comparing  Figure  23  with  Figure  12. 
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In  considering  the  effects  of  the  shock  wave  on  turbulence  which  passes 
through  it,  the  shock  wave  is  generally  assumed  to  be  stationary.  The  factors 
which  are  expected  to  affect  the  turbulence  are,  then,  the  pressure  jump,  density 
jump,  and  velocity  jump  across  the  shock.  However,  it  is  clear  that,  except  at  the 
side  boundaries  where  the  shock  location  is  fixed  by  the  boundary  conditions,  the 
shock  must  move  since  the  shock  speed  is  a  function  of  the  upstream  flow 
conditions.  Since  the  conditions  upstream  of  the  shock  are  varying  from  time-step-to- 
time-step  and  point-to-point,  the  shock  speed  must  vary  locally  as  this  nonuniform 
flow  piisses  through.  In  addition  to  this  fluctuation  of  the  shock  speed,  the  shape 
of  the  shock  must  also  undergo  some  variance  in  accordance  with  the  spatial 
nonuniformity  of  the  flow  field.  This  raises  the  question  of  whether  the  movement 
and  rippling  of  the  shock  should  be  accounted  for  when  calculating  the  statistics  for 
the  turbulence  downstream  of  the  shock  wave. 

Consider  the  equation  for  the  conservation  of  meiss  in  a  compressible  flow 

Pt  +  (/’ujx  +  {P^)y  —  0  (17) 

where  the  subscript  notation  is  used  to  denote  partial  derivatives.  At  a  normal 
shock  wave,  because  of  the  discontinuous  flow  conditions,  the  equation  for  the  jump 
in  conditions  at  the  shock  becomes 

[/)]ft  4-  (H^x  +  (^v]fy  =  0  (18) 

where  the  square  brackets  represent  the  discontinuous  changes  in  the  quantities,  and 
the  shock  location  is  given  by 


f(x,  y,  t)  =  0 

For  steady  flow,  this  reduces  to 


[^u]  =  A  (19) 

where  A  is  the  error  in  the  numerical  solution.  Dividing  Equation  (18)  by  f^  and 
recognizing  that 
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(20) 


fjfx  =  -Us 


and, 


fy/f,  =  -(dXjdy) 


(21) 


the  jump  condition  becomes 


-l/’jus  -h  l/’u]  -  (/)v(dXs/dy)]  =  0  (22) 


Let 

p  -  <p>  +  p' 

u  =  <u>  +  u' 

V  =  <v>  +  v* 

M  =  dX3/dy 

where  <>  represents  the  mean  value  of  each  quantity,  as  previously  defined.  Then 
Equation  22  becomes,  after  averaging  over  the  entire  solution 

-[</JUs>]  +  [<p><u>  +  </)'u'>]  -  1</7vM>]  =  A  (23) 


In  order  to  quantify  the  shock  motion  from  the  results  of  the  numerical 
calculations,  the  shock  location  was  assumed  to  be  the  point  where  the  local 
pressure  passed  through  0.85  with  positive  gradient.  This  was  a  nominal  value 
determined  by  trial  and  error  to  locate  a  point  which  was  consistently  found  in  the 
shock  capture  region. 
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The  change  in  the  turbulence  quantity  <p'\i'>  across  the  shock  is  of  the 
same  order  as  the  changes  in  the  shock  motion  and  shock  ripple  terms  as  shown  in 
Figures  24  and  25.  Both  of  the  shock  terms  undergo  a  strong  perturbation  at  the 
shock  and  then  gradually  decline  as  the  correlation  between  the  turbulence  and  the 
shock  motion  decreases  with  increasing  distance  from  the  shock.  The  reason  for  the 
apparent  negative  correlation  between  the  shock  ripple  and  the  pv  product 
downstream  of  the  shock  is  not  understood  at  this  time.  In  contrast  to  this 
rapidly  declining  correlation,  the  <p'u'>  correlation  declines  much  more  slowly  after 
the  shock  perturbation. 

The  magnitudes  of  the  shock  motion  and  shock  ripple  terms  are  found  to  be 
as  follows: 


<u  ^  0.0556 

NUg  y' 

<U'2^1/2  ^  0.0575 

Both  of  these  quantities  are  roughly  twice  the  magnitude  of  the  turbulence  intensity 
terms  in  the  vicinity  of  the  shock  shown  in  Figure  9. 

It  is  recalled  that  the  shock  capture  region  was  considered  to  be  a  region 
which  could,  be  ignored  with  respect  to  the  numerical  solution,  taking  only  the 
beginning  and  end  points  of  the  region  as  accurate.  Confirmation  of  this 
assessment  is  offered  by  Figure  26  wherein  is  shown  the  net  mass  flux  as  given  by 


<N>  =  -  </3Us>  +  <pu>  -  </5vM> 


(24) 


The  beginning  and  end  of  the  shock  region  are  at  approximately  x  =  0.0275  and 
0.035  respectively,  as  discussed  previously.  The  value  of  the  balance  is  nearly  the 
same  at  the  two  points,  even  though  the  balance  changes  significantly  between  the 
points.  Thus,  the  mass  balance  after  the  shock  is  nearly  the  same  as  it  was  ahead 
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of  the  shock.  This  is  as  it  should  be,  allowing  for  some  entrainment  of  fluid  due 
to  the  fact  that  the  region  for  which  the  statistics  are  calculated  is  not  bounded  by 
solid  walls. 

In  the  interest  of  turbulence  modeling,  the  effects  of  the  shock  wave  on  the 
turbulence  kinetic  energy  is  of  great  importance.  The  quantity  presented  in  Figures 
11  and  22  is  the  kinetic  energy  of  the  turbulence  fluctuations, 

<k>  =  «pu'u’>  +  <pv'v'»/2  (25) 

This  quantity  contains  the  density  and  its  fluctuations  as  well  as  the  velocity 
fluctuations  of  the  turbulent  flow.  It  is  important  to  determine  whether  the 
increase  in  this  quantity  downstream  of  the  shock  observed  in  Figures  11  and  22  is 
due  to  factors  in  addition  to  the  25-percent  increase  in  density  across  the  shock 
jump.  For  this  purpose,  it  is  more  instructive  to  examine  a  quantity  which 
contains  only  turbulence  quantities.  For  turbulence  models,  the  quantity  of  most 
interest  in  the  present  context  would  be 


k*  =  (<u'u'>  +  <v'v'>)/2 


(26) 


The  variation  of  <u'u'>  and  <v'v'>  and  of  k  are  shown  in  Figures  27  and 
28,  respectively.  The  <u'u'>  and  <v'v'>  exhibit  the  same  kind  of  behavior  as 
the  <u'u'>^^^  and  <v'v'>^^^  of  Figure  9.  Similarly,  the  k  variation  indicates  an 
increause  over  the  no-shock  ciise  of  approximately  30-percent.  It  is  also  observed  in 

Figures  9  and  27  and  11  and  28  that  the  decay  rate  of  the  various  quantities  is 

changed  significantly  immediately  downstream  of  the  shock,  increasing  slightly  for 

<v'v'>  and  decreasing  greatly  for  <u'u'>.  This  condition  exists  only  for  a  short 

distance  downstream  of  the  shock  where  the  decay  rate  resumes  a  trend  nearly 
parallel  to  the  no-shock  case.  Thus,  the  kinetic  energy  of  the  turbulence  is 
increased  over  the  no-shock  case  through  an  increase  in  the  intensity  of  the  velocity 
fluctuations  through  the  shock,  and  then  becomes  further  displaced  from  the  no¬ 
shock  case  through  the  apparent  decreEise  in  the  dissipation  and  remains  at  a  higher 
level  than  for  the  no-shock  case  as  the  decay  returns  to  the  unshocked  rate. 
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If  an  analysis  like  that  discussed  previously  with  respect  to  the  conservation  of 
mass  is  performed  for  the  conservation  of  momentum  and  the  equation  for  the 
transport  of  turbulence  kinetic  energy,  it  is  clear  that  terms  containing  the  shock 
motion  and  shock  ripple  will  appear.  It  has  been  demonstrated  that  the  shock 
motion  and  ripple  produce  terms  of  the  same  order  as  the  density-velocity 
correlations  and  turbulence  intensity  of  the  flow  field.  While  an  exact  cause-and- 
effect  relationship  has  not  yet  been  established,  it  is  strongly  suggested  by  these 
observations  that  the  shock  motion  may  be  a  significant  factor  in  the  production  of 
turbulence  kinetic  energy  by  the  shock.  An  analysis  relating  the  increase  in 
turbulence  kinetic  energy  to  the  shock  motion  is  presented  in  Appendix  A. 

Finally,  the  shock  wave/turbulence  interaction  is  observed  to  produce  a  net 
increase  in  the  streamwise  mass  flux  within  the  channel  defined  by  the  turbulent 
inflow.  In  Figure  17,  the  mass  flux  in  that  channel  was  shown  for  the  case  with 
no  shock  wave.  The  mass  flux  increases  slightly  in  the  streamwise  direction. 
Recalling  Figure  16,  where  the  difference  in  the  mass  flux  between  the  shocked  and 
non-shocked  case  is  shown,  the  difference  upstream  of  the  shock  is  essentially  zero, 
except  for  the  oscillations  which  increase  as  the  shock  is  approached  and  then 
decrease  as  the  shock  is  passed.  After  the  shock,  as  the  oscillations  diminish,  the 
difference  shows  an  increasing  trend. 

The  observations  discussed  herein  lead  to  the  conclusion  that  there  is  a 
correlation  between  the  fluctuating  density  in  the  vicinity  of  the  shock  and  the 
fluctuating  local  shock  speed,  and  that  the  simple  jump  relation  that  exists  for  the 
unperturbed  mean  flow,  namely  =  P2^2  does  not  carry  over  to  the  turbulent 

flow.  In  fact,  it  is  found  that  the  correlation  between  p'  and  Ug  is  of  the  same 
order  as  for  p'  and  u‘.  Furthermore,  the  magnitude  of  the  rms  fluctuation  in  Ug 
and  in  dXg/dy  are  both  of  the  same  order  as  the  turbulence  intensity  components. 
A  similar  analysis  may  be  constructed  for  the  momentum  and  energy  equations. 
This  suggests  that  an  important  ingredient  of  a  model  of  turbulence  interacting  with 
a  shock  is  the  coupling  between  the  shock  motion  and  the  density,  energy,  and 
velocity  fields. 
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INDICIAL  ANALYSIS 


In  unsteady  transonic  aerodynamics,  it  is  observed  that  the  shock  motion  is 
nearly  linearly  dependent  on  the  unsteady  forcing  function  and,  because  of  this 
linearity,  an  "indicial"  formulation  of  the  shock  motion  as  proposed  by  Nixon  (Ref. 
5)  can  be  used.  In  such  a  formulation  the  shock  location  for  any  time-dependent 
forcing  function  is  related. to  its  "indicial  response"  to  a  step  change  in  the  forcing 
function  through  a  convolution  integral.  Thus,  if  6Xs(t)  is  the  shock  location,  then 

5X^(t)  -  5X^^(t)€(0)  (27) 

where  (5X3g(t)  is  the  transient  response  to  a  unit  step  change  in  some  parameter  e; 
e(t)  is  an  arbitrary  schedule  of  the  parameter  e,  and  (5Xs(t)  is  then  the  shock 
location  due  to  this  schedule.  In  view  of  the  power  of  this  theory  in  unsteady 
transonic  flow,  it  was  decided  to  investigate  its  applicability  to  turbulence.  That  is, 
the  parameter  e  is  taken  to  be  the  independent  subset  of  the  turbulence  quantities 
u',  v',  /)',  and  p';  at  this  stage  the  correct  subset  has  not  been  determined.  The 
validation  of  this  indicial  function  idea  is  discussed  in  Appendix  B. 

In  the  indicial  analysis  for  transonic  flow  it  is  frequently  convenient  to 
represent  the  numerically  generated  indicial  response  function,  (5Xs£(t),  by  a  simple 
series,  such  as 


61 


(n) 


(t)  = 


(”) 

^  (n) 


[l  -  exp(-a("^t)]  [l 


(28) 


where  (SXgg^"^  is  the  indicial  response  due  to  a  step  change  in  the  variable  of 

6(3^"^  5X5^"^(«')  is  the  value  of  5X3(t)  as  t-**  and  a^"^  and  bj^"^  are  constants 
associated  with  the  variable 
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Consider  now  the  response  of  the  shock  wave  to  the  nth  turbulence  quantity 
f^"^(xj,t)  where  Xj  is  the  spatial  vector  and  t  is  time.  f(xj,t)  may  be  expanded  as  a 
Fourier  series  to  give 


00 

=  Y1  (x.)exp(iwmt) 

■>  m=-“ 


(29) 


where  U  is  some  fundamental  frequency  and 


J  (n) 


is  a  function  of  x. 

J 


If  f^"^(xj,t)  is  identified  with  e  in  Equation  (27),  and  if  the  frequencies  of 
turbulence  are  high,  then  using  the  result  obtained  in  Reference  5,  it  may  be 
shown  that 
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Differentiation  with  respect  to  time  gives  after  some  rearrangement 
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and  thus  the  shock  speed  becomes  a  linear  combination  of  the  turbulence  variables 
f^’'^(Xj,t).  It  can  also  be  seen  from  Equation  (30)  that  if  the  frequency,  u,  of  the 
turbulence  tends  to  infinity,  the  shock  location  does  not  change  but  the  shock  speed 
has  a  magnitude  of  the  same  order  as  the  turbulence  quantities. 


It  is  of  interest  to  consider  the  spatial  derivative  of 
of  Equation  (31)  with  respect  to  xj  gives 
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If  Sf^^^/Qxj  is  the  same  order  of  magnitude  as  Qf^'^^/Qt  then  the  shock  ripple, 
M^“^(xj,t),  has  the  same  order  of  magnitude  as  the  shock  speed,  Ug^"^(xj,t)  and, 
hence,  also  as  the  turbulence  quantities. 

In  this  analysis,  the  constants  a^"^,  and  6X^“^(«>)/e^”^  may  be  dependent 

on  the  geometry  of  the  problem,  a  result  inferred  from  work  in  transonic 
aerodynamics.  Hence,  the  relation  between  the  shock  speed  and  ripple  and  the 
turbulence  quantities  may  not  be  universal. 

The  analysis  in  this  section  implies  the  following  relations. 

(a) .  The  shock  speed  is  the  same  order  of  magnitude  as  the  turbulence 

quantities  and  is  a  linear  combination  of  them. 

(b) .  The  shock  ripple  may  be  of  the  same  order  of  magnitude  as  the  shock 

speed. 

(c) .  The  shock  motion  is  an  order  of  magnitude  smaller  than  the  shock 

speed. 
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CONCLUSIONS  AND  RECOMMENDATIONS 


The  work  reported  herein  represents  the  first  step  in  a  study  of  shock  waves 
interacting  with  turbulence.  The  overall  goals  of  the  study  are  to  determine  the 
mutual  effects  of  the  interaction  and  to  determine  improvements  for  turbulence 
models  in  the  presence  of  shock  waves.  The  ultimate  goal  is  to  extend  this 
knowledge  to  improving  turbulence  models  in  the  presence  of  shock  waves  at 
hypersonic  speeds. 

Two  short  term  goals  for  the  present  study  were  to  develop  an  accurate 
numerical  approach  in  two  dimensions  which  would  serve  to  illuminate  the  problems 
that  may  be  encountered  in  the  long  range  study  and  to  aid  in  determining 
methods  of  analyzing  the  turbulence.  These  goals  were  accomplished  by  developing 
a  computer  code  employing  MacCormack's  explicit  numerical  method  to  solve  the 
Navier-Stokes  equations  for  a  normal  shock  wave  encountering  turbulence. 

Another  short  term  goal  was  to  develop  insight  into  shock/turbulence 
interaction.  This  was  accomplished  by  employing  a  low  supersonic  Mach  number 
which  allowed  the  study  to  take  place  in  an  environment  for  which  the  basic 
elements  of  the  flow  are  well  understood,  and  also  to  make  use  of  certain 
expedients  available  from  transonic  flow  theory.  The  major  insight  that  was 
obtained  from  the  work  was  that  shock/turbulence  interaction  contains  influences 
from  the  shock  motion  that  are  quite  unapparent  from  a  Reynolds  averaged  point 
of  view.  This  conclusion  must  be  examined  further  using  three-dimensional 
turbulence  and  higher  Mach  numbers.  Also,  the  work  must  be  extended  to 
encompass  oblique  shocks  interacting  with  turbulence  in  a  shear  flow,  since  such  are 
the  flows  that  will  ultimately  be  of  most  interest. 

Another  short  term  goal  of  the  work  was  to  extend  the  numerical  approach  to 
three  dimensions.  The  effort  to  date  has  centered  on  developing  the  appropriate 
analytical  and  computational  tools  for  the  analysis  of  turbulence/shock  wave 
interactions.  The  computer  code  used  in  these  preliminary  calculations  has  been 
combined  with  another  code  employing  the  same  algorithm  already  possessing  the 
required  three-dimensional  capability.  The  turbulence  statistics,  boundary  condition 
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treatment,  and  turbulence  generation  used  in  the  two-dimensional  code  have  thus 
been  extended  to  three  dimensions.  The  code  is  now  ready  for  a  fully  three- 
dimensional  treatment  of  the  turbulence. 

Examination  of  the  solutions  generated  by  the  two-dimensional  code  indicated 
that  the  numerical  errors  that  are  inherent  in  the  explicit  MacCormack  algorithm 
do  not  affect  the  turbulence  calculations.  Large  oscillations  are  produced  by 
discontinuous  changes  in  the  mean  flow  quantities,  but  the  perturbation  quantities 
undergo  smaller  changes  and  change  more  smoothly.  Therefore,  the  pertubation 
quantities  exhibit  only  relatively  small  oscillations.  Thus,  the  turbulence  calculations 
are  believed  to  be  accurate  to  within  the  truncation  error  for  the  numerical  method. 

Several  observations  were  made  which  indicate  the  shock  has  a  significant 
effect  on  the  turbulence.  The  shock  produces  a  jump  in  the  turbulence  quantities, 
with  a  long  relaxation  distance  to  return  to  unshocked  values.  The  turbulence 
kinetic  energy  is  increased  by  30  percent  by  the  shock.  The  density-velocity 
correlation,  </3'u'>,  becomes  important  during  the  shock  jump  and  is  greatly 
increased  over  the  case  without  a  shock.  On  the  other  hand,  the  pressure-velocity 
correlation,  <p'u'>,  is  not  quite  so  important. 

Using  the  no-shock  case  as  a  guide,  the  rate  of  decay  before  and  after  the 
shock  indicates  that  the  dissipation,  e,  changes  only  slightly  compared  to  the  kinetic 
energy.  Thus,  if  eddy  viscosity  is  defined  as  i't  -  k  /e,  the  ratio  of  i/^  before  and 
after  the  shock  is  kj  /k2  -  2.  This  has  implications  regarding  heat  transfer 

which  can  change  dramatically  because  of  changes  in 

The  results  show  that  shock  speed  and  ripple  may  be  important  factors  in 
determining  the  turbulence  downstream  of  a  shock  wave.  Shock  speed  and  ripple 
correlations  are  the  same  size  as  other  important  turbulence  correlations,  such  as, 
<p'u'>,  <u'u'>,  <p'u'>,  and  k;  shock  ripple  terms  may  be  more  important  for 
oblique  shocks  because  of  larger  v  component. 

Transonic  indicial  theory  applied  to  the  shock/turbulence  interaction  leads  to 
several  interesting  observations.  For  example,  while  undergoing  fluctuations  in 
velocity  equal  to  the  turbulent  velocity  fluctuations  of  the  flow  field,  the  shock 
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location  moves  only  slightly.  Another  observation  is  that  shock  speed  correlations 
are  a  linear  combination  of  correlations  of  turbulent  flow  quantities.  Jump 
correlations  for  shock  speed  are  about  the  same  order  of  magnitude  as  the  flow 
correlation  jumps. 

In  conclusion,  this  work  has  demonstrated  that  there  are  potentially  very 
important  interactions  between  a  shock  wave  and  a  turbulent  flow  field.  The 
approach  taken,  employing  available  numerical  computational  tools  and  an 
approximate  representation  of  turbulence,  has  provided  useful  insights  for  expanding 
the  research  to  encompass  more  realistic  flow  configurations  and  conditions.  While 
two-dimensional  turbulence  is  not  realistic,  for  the  purposes  outlined  in  this  report 
it  contains  the  essential  elements  of  a  fluctuating  flow  field.  It  is  strongly 
recommended  that  future  calculations  be  done  in  a  three-dimensional  turbulence 
field.  Also,  a  larger  computational  domain  is  needed  to  check  the  relaxation 
distiuice  of  the  shock  effects  on  the  turbulence  and  to  isolate  the  computation  from 
boundary  effects.  A  more  accurate  algorithm  is  needed  in  order  to  eliminate  any 
question  of  the  effect  of  the  numerical  errors  on  the  observed  turbulence  after  the 
shock.  Finally,  calculations  should  be  done  for  an  oblique  shock  in  order  to  obtain 
stronger  v  correlations. 
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Table  I. 


List  of 

rms  values  of  u 

'  .  v'  ,  p' ,  and 

p'  for  steady 

shock  case . 

X 

u '  rms 

V  '  rms 

p '  rms 

p  '  rms 

0 . 00000 

7 . 12319E-06 

1 . 44329E-06 

7 . 39986E-06 

7.03310E-06 

0.00125 

2 . 44887E-06 

0 . OOOOOE+00 

5 . 96046E-08 

1 .04606E-06 

0 . 00250 

7 . 20255E-06 

1.44321E-06 

7 . 37847E-06 

7 . 02020E-06 

0 . 00375 

3.31865E-06 

7 . 10731E-08 

3 . 00102E-06 

2.51402E-06 

0 . 00500 

8 . 46386E-06 

2 . 27635E-06 

8.73325E-06 

8.37122E-06 

0.00625 

5 . 46286E-06 

4 . 38622B-07 

5 . 15468E-06 

5.04989E-06 

0 . 00750 

1.04541E-05 

3.45727E-06 

1 .07116E-05 

1 .03719E-05 

0.00875 

8 . 53908E-06 

1 .38007B-06 

8.41271E-06 

8.31928E-06 

0.01000 

1.31752E-05 

5 . 23048E-06 

1.34059E-05 

1.31232B-05 

0.01125 

1 . 26839E-05 

3.21078B-06 

1 . 27408E-05 

1 .26148E-05 

0.01250 

1.75992E-05 

7.98542E-06 

1 .79207E-05 

1 . 75934E-05 

0.01375 

1 .91810E-05 

6 . 38178B-06 

1.95491E-05 

1.93972E-05 

0.01500 

2.51577E-05 

1 . 20906B-05 

2 . 58342E-05 

2.64843E-05 

0.01625 

2.88457E-05 

1 . 10S56E-05 

2 . 97493E-05 

2 . 95925E-05 

0.01750 

3 . 50163E-05 

1.69789E-05 

3 . 62608E-05 

3.59802E-05 

0.01875 

3.80046E-05 

1.54S31E-05 

3 . 94589E-05 

3.94713E-05 

0.02000 

4.06055E-05 

1.92527E-05 

4.20449E-05 

4.20517E-05 

0.02125 

3.80670E-05 

1.47329E-05 

3.92283E-05 

3.95789B-05 

0.02250 

3.49481E-05 

1.47206B-05 

3.55554E-05 

3 . 57724E-05 

0.02375 

2 . 90470E-05 

7.81S80E-06 

2 . 92946B-05 

2.96241B-05 

0.02500 

3.08264E-05 

8.70620B-06 

3.08959E-05 

3.08293E-05 

0.02625 

2.90521E-05 

4.04505E-06 

2.96802E-05 

2.95422E-05 

0.02750 

5 . 16982E-05 

8.59280E-06 

5.07960E-05 

5 . 16102E-05 

0.02875 

7.75456E-06 

1.03651E-06 

1 . 24250E-05 

9.88898E-06 

0.03000 

2 . 23273E-04 

1 . 13871E-05 

2.26297E-04 

2.37557E-04 

0.03125 

1 .61209E-04 

4.52128E-05 

1 . 95533E-04 

2 . 04737E-04 

0 . 03250 

3.78583E-05 

3.66188E-05 

3.83401E-05 

4.57391E-05 

0 . 03375 

1 .80216E-05 

3.88444E-05 

1.92622E-05 

2.06197E-05 

O . 03500 

5.37351E-05 

2.77664E-05 

6 . 0g432E-05 

6 . 77737E-05 

0 . 03625 

2 . 63460E-05 

2 . 20543B-05 

2 . 66304E-05 

2 . 92868B-05 

0.03750 

4 . 46592E-05 

1.36363E-05 

5 . 04349E-05 

5.59164E-05 

0 . 03875 

2.65819E-05 

7.22448E-06 

2.80054E-05 

3.07446E-05 

O . 04000 

3 . 18990B-05 

3 . 29671E-06 

3 . 63368E-05 

4.03315E-05 

0.04125 

2 . 05268B-05 

6 . 16266E-06 

2 . 25863E-05 

2.47881E-05 

0.04250 

1 .g4g3gE-05 

7.88442E-06 

2.26921E-05 

2.53252E-05 

0 . 04375 

1 . 57339E-05 

1 . 13655E-05 

1 . 90377E-05 

2.0gi83E-05 

0 . 04500 

1 . 16217E-05 

1 . 17g25E-05 

1 .37166E-05 

1 .53833E-05 

0 . 04625 

1 . 52566E-05 

1.35844E-05 

1 .g8255E-05 

2 . 18449E-05 

0 . 04750 

1 .01664E-05 

1 . 30935E-05 

1 . 12954E-05 

1 .24460E-05 

0 . 04875 

1 .47201E-05 

1.37006E-05 

1 ,g5935E-05 

2. 16414E-05 
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Table  II. 


List  of  <u>,  <v>,  <^>,  <p>,  and  </5u>  for  steady  shock  case. 


Theoretical 

Values  at  shock 

jump: 

“1 

=  1.15 

U2 

=  0.916305 

PI 

=  1.0 

P2 

=  1.25504 

^1 

=  0.0 

V2 

=  0.0 

PI 

=  0.714286 

P2 

=  0.98286 

<fiu>i 

=  1.15 

<fiu>2 

=  1.15 

X 

<u> 

<v> 

<P> 

<P> 

</>u> 

0.00000 

1.150027 

6.07867E-08 

0.999971 

0.714258 

1.14999 

0.00125 

1.150000 

O.OOOOOE+00 

1.000000 

0.714286 

1 . 15000 

0.00250 

1.149972 

-6.07837E-08 

1.000029 

0.714314 

1.15001 

0.00375 

1 . 150010 

1.02647E-11 

0.999991 

0.714276 

1.15000 

0.00500 

1 . 149964 

-9.57456E-08 

1.000037 

0.714322 

1.15001 

0.00625 

1 . 150024 

1.39344E-08 

0.999977 

0.714262 

1 . 15000 

0.00750 

1 . 149947 

-1.45508E-07 

1.000054 

0.714339 

1 . 15001 

0.00875 

1.150047 

5.19740E-08 

0.999953 

0.714238 

1.14999 

0.01000 

1.149919 

-2,23606E-07 

1.000083 

0.714367 

1.15001 

0.01125 

1.150086 

1.35955E-07 

0.999913 

0.714199 

1.14999 

0.01250 

1.149871 

-3.62275E-07 

1.000132 

0.714415 

1 . 15002 

0.01375 

1.150150 

3.14396E-07 

0.999848 

0.714134 

1 . 14999 

0.01500 

1.149790 

-6.27321E-07 

1.000214 

0.714497 

1 . 15004 

0.01625 

1 . 150258 

6.64969E-07 

0.999738 

0.714026 

1.14996 

0.01750 

1.149652 

-1.07788E-06 

1.000355 

0.714636 

1 . 15006 

0.01875 

1.150439 

1.16248E-06 

0.999554 

0.713844 

1.14993 

0.02000 

1 . 149418 

-1.51634E-06 

1.000592 

0.714871 

1.15010 

0.02125 

1 . 150735 

1.36886E-06 

0.999252 

0.713545 

1 . 14987 

0.02250 

1 . 149020 

-1.35362E-06 

1.000995 

0.715272 

1.15016 

0.02375 

1.151165 

8.29798E-07 

0.998808 

0.713109 

1 . 14979 

0.02500 

1.148177 

-6.89811E-07 

1.001815 

0.716110 

1 . 15026 

0.02625 

1.151036 

3.19065E-07 

0.998812 

0.713187 

1.14967 

0.02750 

1 . 143497 

-5.40899E-07 

1.006029 

0.720638 

1.15039 

0.02875 

1.137744 

1.03004E-07 

1.010291 

0.725830 

1 . 14945 

0.03000 

1.084217 

-7.05174E-07 

1.061083 

0.780058 

1 . 15044 

0.03125 

0.961642 

-2.69602E-06 

1.195909 

0.922350 

1 . 15004 

0.03250 

0.924830 

-2.22471E-06 

1.243277 

0.971103 

1.14982 

0.03375 

0.918577 

-2.15126E-06 

1.251994 

0.979883 

1 . 15005 

0.03500 

0.916195 

-1.53559E-06 

1.255089 

0.983130 

1.14991 

0.03625 

0.916602 

-9.97246E-07 

1.254663 

0.982625 

1 . 15003 

0.03750 

0.915980 

-5.79323E-07 

1.255409 

0.983443 

1.14993 

0.03875 

0.916406 

-5.55969E-08 

1.254908 

0.982887 

1 . 15000 

0.04000 

0.916056 

1.42036E-07 

1.255322 

0.983344 

1.14995 

0.04125 

0.916332 

4.97701E-07 

1.254997 

0.982985 

1.14999 

0.04250 

0.916121 

5.41513E-07 

1.255247 

0.983261 

1.14996 

0.04375 

0.916295 

7.36329E-07 

1.255042 

0.983035 

1.14999 

0.04500 

0.916169 

6.99422E-07 

1.255191 

0.983200 

1.14997 

0.04625 

0.916278 

7.83975E-07 

1.255063 

0.983058 

1.14999 

0.04750 

0.916205 

7.14981E-07 

1.255150 

0.983155 

1 . 14998 

0.04875 

0.916273 

7.35186E-07 

1.255071 

0.983067 

1.14999 
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Figure  2.  Convergence  history  of  inaxiinum  error 


-■{6- 


10000;  Second  order  artificial  dissipation. 


Figure  3.  Continued 
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c.  Re  =  10000;  Fourth  order  artificial  dissipation. 


Figure  3-  Continued 


d.  Re  =  10000;  Second  and  fourth  order  artificial  dissipation. 


!•  ■..c.iirc  3.  Continued 


Re 


100000,  Second  and  fourth  order  artificial  dissipation. 


Figure  3  Continued 


1000;  Second  and  fourth  order  artificial  dissipation. 


Figure  3.  Concluded 
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Figure  9.  Statistics  of  turbulent  velocity  field  through  a  shock. 
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Figure  10.  Statistics  of  pressure  and  density. 
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Figure  15.  Plot  of  <pu>  for  turbulent  flow  with  a  shock  wave. 


APPENDIX  A 


ANALYSIS  OF  TURBULENT  KINETIC  ENERGY 
Consider  the  conservation  of  mass  through  a  normal  shock  wave,  that  is 

where  p  is  the  density  and  u  is  the  velocity  normal  to  the  shock;  the  subscripts  1 
and  2  denote  values  ahead  of  and  behind  the  shock  respectively. 

The  density  can  be  written  in  the  form 

P  =  ^  - 2 -  «p(-s/R)  (A2) 

where  u  is  the  velocity  normalized  with  respect  to  a  reference  velocity,  uo6,  with  an 
associated  Mach  number,  M®;  S  is  the  entropy  and  R  is  the  gas  constant. 

If  the  shock  wave  is  weak,  the  entropy  jump,  AS,  across  the  shock  is  given  by 

27  2  1  2  3 

AS/R  =  - r  (M/  -  =  A(M/  -  1)'*  (A3) 

3(7-1)"  ^  ^ 

where  Mj  is  the  Mach  number  ahead  of  the  shock  wave. 

Now  let  u  be  represented  by 

u  =  1  +  u  (A4) 

Using  Equations  (A2),  (A3),  and  (A4),  Equation  (1)  can  be  expanded  to  second 
order  in  u  as 

1  .  p\  -  f  =  (1  .  p\  -  I  u/)  [1  -  A(Mj2  _  ^)3^  (^5^ 

where 


A1  - 


(A6) 


=  1  - 

=  m/[3  +  (7-2)m/] 

In  this  formulation  Mj^  is  given  by 

Mi^  =  (A7) 

The  perturbation  velocity  can  now  be  represented  by 

u  =  u  +  u'  (A8) 

where  u  is  a  stationary  value  and  u '  is  a  fluctuating  value. 

Now  let 


u^  =  0 

and  a  zeroth  order  perturbation  analysis  gives 

fl2-  (C  -  2  ,fl6  n 

P  U2  -  2  ^2  *  ® 

while  a  first  order  analysis  gives 

=  (P%'  -  «U2’U2')  {1  +  [-P®  -  + 

If  the  shock  is  weak  such  that 

<  <  1 

then  a  first  approximation  to  Equation  (All)  is 

P^uj'  -  =  P^ug'  -  «U2'U2' 


(A9) 


(AlO) 


(All) 


(A12) 


(A13) 


A2  - 


Now  the  above  analysis  assumes  that  the  velocities  are  relative  to  the  shock  wave. 
If  the  shock  wave  is  moving  at  speed  u'j  and  is  at  an  inclination  m,  where 


then,  in  a  fixed  frame  of  reference 


u‘  =  (u*  -  u  ‘  -  mv  -  +  m  ) 


2,1/2 


(A15) 


where  u  is  the  fluctuating  velocity  in  the  fixed  frame  of  reference,  v  is  the  mean 
value  of  V  and  v'  is  the  fluctuating  value  of  v.  Using  Equation  (A15)  in  Equation 
(A  13)  3md  taking  averages 

x[<U2^'u2'>  -  <Ug'u2'>]  +  =  /c[<U2'u2'>  -  /?^<mv2'>  (A16) 

where  terms  of  order  higher  than  the  second  have  been  neglected.  Equation  (A16) 
can  be  rearranged  to  give 


<uj'u2'>  =  <U2'u2'>  +  [<Ug'uj'>  -  <Ug'u2’>]  -  P  1^  [<mV]^'>  -  <mv2'>] 

If  the  shock  is  stationary,  then  Equation  (A17)  reduces  to 


(A17) 


<u,  'u,  •>  =  <U,,'Uo'> 


(A18) 


For  the  tangential  component  of  velocity  the  jump  in  v  in  shock  fixed 
coordinates,  gives 


'"l  ■  ^2 


(A19) 


and  in  space  fixed  coordinates 


Vi  +  (uj^  -  Ug)m  =  V2  +  (U2  -  Ug)m 

Again  splitting  v  into  a  stationary  and  non-stationary  component  gives 


(A20) 


A3  - 


Vj'  +  (uj'  -  Ug')m  =  V2'  +  (U2'  -  Ug)iii  (A21) 

squaring  both  sides,  neglecting  terms  of  order  higher  tlian  the  second  and  averaging 
gives 

<V2'vj'>  =  <V2‘v2'>  (A22) 

Using  Equation  (A17)  and  (A22)  it  may  be  seen  that  the  jump  in  turbulent  kinetic 
energy,  defined  by 


k  =  J  {<u'u'>  +  <v'v'>}  (A23) 

is  given  by 

~  ^^2  ^  2  { "  <mv2'>]}  (A24) 

and  it  may  be  noted  that  if  the  shock  speed  and  ripple  terms  are  neglected,  then  k 
is  continuous  through  a  shock  wave. 

If  the  data  from  the  simulations  are  used  to  evaluate  the  right  hand  side  of 
Equation  (A24)  then 


ki  -  k2  =  -2.31  X  10  (A25) 

which  is  an  increase  of  turbulent  kinetic  energy  through  the  shock.  In  the 
simulation,  k,  as  defined  in  Equation  (23),  is  5.56  x  10'^  and,  hence,  the  shock 
increases  k  by  41%. 
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APPENDIX  B 

INDICIAL  ANALYSIS  ' 

i 

Since  the  flow  conditions  being  used  in  these  investigations  are  in  the 
transonic  range,  it  was  suggested  that  some  insight  for  assessing  the  magnitude  of 
the  effect  of  perturbations  on  the  shock  speed  and  the  influence  of  the  shock  speed 
on  the  turbulence  might  be  possible  from  a  linear  theory  such  as  the  transonic 
indicial  theory  developed  by  Nixon  (Refs.  Bl  and  B2). 

The  indicial  method  uses  the  principle  of  superposition  to  relate  complex 
motions  to  an  indicial  response  by  means  of  Duhamel's  integral.  The  indicial 
response  is  the  time  history  of  the  motion  in  question  in  response  to  a  step  change 
in  some  flow  parameter. 

The  basic  equation  for  shock  movement  is,  from  Reference  Bl 

6X(t)  =  <yx  (t)e(O)  +  fdX  (r)  dr  (Bl) 

s  Sg  qJ  Sg  on 

where  ffX,  is  the  shock  movement,  and  €(t)  is  a  time-dependent  parameter  typifying 
the  unsteady  motion.  5X3e(t)  is  the  indicial  response,  defined  as  the  response  of  the 
shock  location  due  to  the  instantaneous  unit  change  in  the  parameter  e(t).  In 
Reference  B2,  Nixon  shows  that  for  an  input  representing  simple  harmonic  motion, 
the  shock  motion  is  also  harmonic,  with  a  phase  lag  varying  with  the  frequency  of 
the  driving  function.  In  particular,  for  very  high  frequencies,  the  shock  motion  goes 
to  zero  as  1/k  where  k  is  the  reduced  frequency,  wc/uo#,  where  c  is  a  reference 
length  and  u«  is  the  freestream  velocity.  However,  while  the  magnitude  of  the 
shock  movement  vanishes  as  frequency  increases,  the  theory  indicates  that  the 
velocity  of  the  motion  remains  oscillatory  with  a  finite  amplitude  while  the 
amplitude  of  the  acceleration  of  the  shock  motion  approaches  «  as  the  reduced 
frequency. 

The  above  discussion  suggests  that  an  important  part  of  the  turbulence  in  a 
shock/turbulence  interaction  may  be  the  motion  of  the  shock  itself,  which  is  not 
accounted  for  by  the  Reynolds  averaging  process  which  averages  the  fluctuating 
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quantities  over  a  long  time  and  considers  only  fluctuations  of  the  flow  velocities  and 
other  flow  variables  relative  to  a  fixed  coordinate  system.  In  order  to  obtain  an 
estimate  of  the  magnitude  of  the  shock  motion  corresponding  to  the  numerical 
experiments  discussed  previously,  a  series  of  calculations  were  made  for  a  case  using 
the  indicial  theory. 

In  the  indicial  test  cases,  the  same  type  of  flow  configuration  was  used  as  for 
the  results  shown  in  Figures  9-17.  However,  instead  of  imposing  a  random 
variation  of  the  flow  variables  in  the  center  portion  of  the  inflow  boundary,  two 
cases  were  calculated.  The  first  case  was  the  indicial  response  case  in  which  all 
flow  variables  were  held  constant  at  their  initial  values  on  the  inflow  boundary 
except  for  the  streamwise  velocity  component,  u,  which  was  given  a  step 
perturbation  at  time  zero  in  the  center  portion  of  the  boundary  (40  grid  points) 
and  held  at  the  perturbed  value  for  the  entire  calculation.  In  the  second  case,  the 
u-velocity  component  was  calculated  from  the  random  function  in  time,  the  same 
value  being  specified  at  all  y  locations  on  the  inflow  boundary  for  each  time  step. 
In  other  words,  the  inflow  was  a  two-dimensional  pulse  with  random  fluctuations. 

Typical  results  from  the  indicial  response  calculation  are  shown  in  Figures  BI¬ 
BS.  In  Figure  Bl,  contours  of  constant  pressure  are  shown  1000  time  steps  after 
\the  initial  step  perturbation  of  the  inflow  velocity.  In  Figure  B2,  a  composite  plot 
is  shown  of  the  pressure  distributions  along  constant  x  lines  near  the  centerline  of 
the  flow.  The  initial  pressure  distribution  of  the  unperturbed  flow  is  shown  for 
comparison.  Finally,  in  Figure  B3,  the  indicial  shock  displacement  from  the  initial 
location  is  shown.  There  is  a  short  delay  from  the  start  of  the  calculation, 
corresponding  to  the  propagation  of  the  step  pulse  through  the  upstream  flow,  then 
the  shock  begins  to  move  in  response  to  the  new  inflow.  The  waviness  of  the 
shock  motion  is  a  numerical  phenomenon  related  to  the  size  of  the  computational 
grid,  the  small  number  of  grid  points  in  the  shock  capture  region,  and  the  linear 
interpolation  of  pressure  which  was  used  to  locate  the  shock.  For  the  calculations 
shown  here,  the  shock  location  was  taken  to  be  the  point  where  the  pressure  was 
equal  to  0.85  and  was  rising  with  increasing  x. 


The  fluctuating  inflow  velocity  is  shown  in  Figure  B4.  This  function  was 
imposed  at  the  same  40  points  in  the  center  portion  of  the  inflow  boundary  as  was 
the  step  pulse  in  the  indicial  calculation.  Two  typical  pressure  contour  plots  are 
shown  in  Figure  B5  for  this  case.  The  first,  Figure  B5a,  shows  the  situation  after 
200  time  steps,  the  second  shows  the  situation  after  1000  steps.  Finally,  the 
distribution  of  pressure  along  x  lines  near  the  centerline  for  the  perturbed  flow  is 
shown  in  Figure  B6a  and  b  for  the  same  two  time  steps  as  the  previous  figure. 

The  motion  of  the  shock  at  the  centerline  was  determined  for  the  perturbed 
flow  in  the  same  manner  as  for  the  indicial  flow,  and  the  location  of  the  shock  as 
a  function  of  time  is  shown  in  Figure  B7.  It  is  noted  that  the  shock  motion  does 
not  oscillate  about  zero  as  a  mean  value.  The  reason  for  this  is  that  the  input 
velocity  perturbation  (Fig.  B4)  does  not  oscillate  about  zero  as  a  mean.  While  the 
perturbation  of  the  velocity  should  eventually  have  zero  mean  value,  the  nonzero 
value  apparent  in  Figure  B5  is  believed  to  be  due  to  the  small  sample  size.  If  the 
calculation  were  to  be  extended  to  include  10000  time  steps,  the  true  statistical 
mean  should  be  apparent.  The  present  analysis  makes  use  of  the  small  sample  size 
because  the  interest  is  in  the  transient  phenomena,  not  in  the  asymptotic  statistics. 

The  shock  perturbation  velocity  calculated  by  differencing  the  Xg  data  is 
shown  in  Figure  B8.  It  is  noted  that  the  magnitude  of  the  shock  velocity  is  of  the 
same  order  as  that  of  the  inflow  velocity  fluctuations.  The  shock  motion  appears 
smoother  than  the  inflow  velocity,  due  to  the  damping  effects  of  the  viscosity,  both 
real  and  artificial,  which  reduces  the  higher  frequency  fluctuations  in  the  calculated 
flow  field. 

With  the  apparent  conceptual  agreement  between  the  full  Navier-Stokes 
solution  and  the  analytical  indicial  theory  regarding  the  magnitude  of  the  shock 
velocity  and  its  potential  importance  to  the  shock/turbulence  interaction,  an 
examination  was  made  of  the  ability  of  the  indicial  theory  to  predict  the  shock 
motion.  Recalling  Equation  (Bl),  an  equation  for  the  shock  velocity  can  be 
obtained  by  differentiation.  Thus, 
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The  indicial  response  was  shown  in  Figure  B3.  The  function  e(t)  was  taken 
to  be  the  inflow  perturbation  velocity,  u(t),  shown  in  Figure  B4.  Calculation  of  the 
integral  in  Equation  (B2)  by  a  simple  trapezoidal  rule  and  the  derivative  of  the 
indicial  response  by  backward  differences  produced  the  distribution  of  the  shock 
velocity  shown  in  Figure  B9.  Comparison  of  Figures  B9  and  B8  reveals  that  the 
shock  velocity  predicted  by  the  indicial  formulation  has  generally  the  same  frequency 
and  phase  as  that  given  by  the  Navier-Stokes  solution,  with  large  errors  in  the 
amplitude  of  the  velocity  fluctuations;  the  reason  for  the  errors  is  not  clear  at 
present.  The  shock  location  determined  by  integration  of  Figure  B9  compares  fairly 
well  with  the  direct  solution  (Fig.  BlO). 

From  the  surprisingly  good  agreement  between  the  indicial  theory  and  the 
Navier-Stokes  solution  for  the  shock  location,  it  can  be  tentatively  concluded  that 
the  assumption  of  a  linear  relationship  between  certain  flow  properties  and  small 
perturbations  of  velocity,  pressure,  or  density  can  be  a  useful  tool  for  providing 
some  insight  to  even  the  complex  interaction  between  a  shock  and  a  turbulent  flow 
field.  For  the  present,  the  primary  conclusion  that  can  be  drawn  is  that  the  shock 
motion  must  be  included  in  any  analysis  of  turbulence  models  for  flows  with  shocks. 
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Figure  Bl.  Contours  of  constant  pressure  for  indicial  solution. 


Figure  B2.  Prosstirt-  disiributioii  near  centerline  for  indicial  solution. 


Figure  B4.  Plot  of  u  from  random  function. 
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a.  t  =  200. 


Figure  B6. 


Pressure  distribution  near  centerline  for  perturbed  solution. 


b.  t  =  1000. 


Figure  B6.  Concluded 
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